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Abstract 

A multicanonical formalism is introduced to describe statistical equilibrium of complex systems 
exhibiting a hierarchy of time and length scales, where the hierarchical structure is described as a set 
of nested "internal heat reservoirs" with fluctuating "temperatures." The probability distribution 
of states at small scales is written as an appropriate averaging of the large-scale distribution (the 
Boltzmann-Gibbs distribution) over these effective internal degrees of freedom. For a large class of 
systems the multicanonical distribution is given explicitly in terms of generalized hyper geometric 
functions. As a concrete example, it is shown that generalized hypergeometric distributions describe 
remarkably well the statistics of acceleration measurements in Lagrangian turbulence. 

PACS numbers: 05.40.-a, 05.10.Gg, 05.70.Ce, 47.27.eb 
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In this Letter we introduce a general formalism to describe statistical equilibrium of 
complex systems with multiple time and length scales. We adopt a viewpoint akin to the 
canonical ensemble perspective — only slightly augmented — , in which the system is given 
a certain temperature by being placed in an infinite heat bath of the proper temperature 
The specific question we wish to address here is how a small subsystem within the 
main system comes into thermal equilibrium with the heat bath and what is the resulting 
orobability distribution of states for such subsystem. In the standard canonical treatment 
l| , the small subsystem and its large surroundings are assumed to be independent and thus 
described by the same distribution law — the Boltzmann-Gibbs (BG) distribution. There 

n 

are, however, many physical systems, such as highly turbulent flows |2J, where owing to 
the existence of a hierarchy of dynamical structures the relevant distributions depend on 
the scale at which the measurements are made. In such cases, the canonical hypothesis 
must be modified accordingly to take into account the more complex process of energy ex- 
change between the subsystem and the heat bath, which will be mediated by the intervening 
hierarchical structure. 

Here we give an effective description of the dynamical hierarchy in terms of a set of 
nested "internal heat reservoirs," where the innermost reservoir surrounds the subsystem of 
interest while the outermost one is in contact with the external heat bath. The complex 
(intermittent) energy flow between adjacent hierarchical levels is then modelled by allowing 
the "temperatures" of such internal reservoirs to fluctuate according to a stochastic dynam- 
ics described by a deterministic term, given by the usual Newton's law of cooling, plus a 
multiplicative noise. (Without the stochastic term the system would, of course, relax to the 
usual Gibbsian equilibrium.) In such scenario, it turns out that for a large class of systems 
the equilibrium distribution can be written explicitly in terms of certain generalized hyperge- 
ometric functions. This family of generalized hypergeometric (GHG) distributions includes, 
as its first two members, the BG distribution and the g-exponential distribution, also known 

as Tsallis distribution, which has been much studied in the context of the so-called nonex- 

P 

tensive thermodynamics [3] . It is thus argued that GHG distributions of higher order should 
naturally appear in complex systems having more than two distinct time scales. As a con- 



crete example, we show that the GHG distribution of seventh order descri 



well the statistics of acceleration measurements in Lagrangian turbulence {4]. 

We consider a system of size L in contact with an external heat reservoir at a fixed 



Des remarkably 



temperature T . We assume that the system possesses a hierarchy of dynamical structures 
of characteristic sizes £j, where for definiteness we take £j = L/2\ for j = 0, 1, ...,n. It is 
also supposed that there is a wide separation of time scales, with smaller structures having 
faster dynamics. We focus our attention on a small subsystem of size r < l n . (One can think 
of this small subsystem as the measurement volume.) In the standard canonical formalism, 
the large subsystem surrounding the small subsystem can be viewed either as a heat bath or 
as a large collection (ensemble) of small subsystems essentially identical with the subsystem 
in focus. Owing to the presence of multiple scales, neither one of these two viewpoints is 
however applicable in our case. We shall instead regard the large subsystem as consisting 
of a set of nested "internal heat reservoirs," where each such reservoir is characterized by 
its own effective "temperature" Tj, with j = 0,1, ...,n. The temperature Tj represents a 
measure of the average energy (at a given time) in the structures of characteristic size £j and 
as such will be treated as a fluctuating quantity, whose probability density function (PDF) 
will be denoted by f(Tj). 

Our aim here is to obtain the probability, P r (ei), of finding a small subsystem of size r 
in a given state of energy £j. By assumption, the subsystem has a dynamics much faster 
than that of the temperature T n of its immediate surroundings. It is therefore reasonable 
to suppose that before T n changes appreciably the subsystem will reach a quasi-equilibrium 
state described by the BG distribution at temperature T n : 

where k is Boltzmann's constant and 

Z (T) = E ex P = [ 9(E)exp (-^j dE. (2) 

Here the energy E is regarded as a continuous variable and g(E) denotes the density of 
states. The marginal distribution P r (ei) is then given by 

Notice that at the largest scale (i.e., n = 0), Eq. ([3]) reduces to the canonical distribution, 
since the external heat bath is assumed to have a constant temperature: f(T) = 5(T — To). 
The distribution P r (e) above generalizes the canonical distribution for systems with multiple 
scales, and hence it will be called multicanonical j^]. 
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The idea expressed in Eq. d3J) of writing the distribution law at small scales as a mixture 
of large-scale distributions has been extensively used in turbulence, where the large-scale 
distribution is usually a Gaussian and various weighting distributions have been considered 

6- 

perstatistics 



loj |. More recently, this idea has also been applied in the context of the so-called su- 
I2I of nonequilibrium systems, where one considers the superposition of 
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two statistics describing the stationary state and a slowly fluctuating intensive parameter, 
usually the inverse temperature (3 = 1/kT. Other approaches related to superstatics have 
also been considered [13L Il4 ]. Although the multicanonical distribution ([3]) bears a formal 
resemblance to superstatistics, our formalism is conceptually quite distinct. In particular, 
we do not prescribe a priori the weighting distribution f(T n ) but rather derive it from a 
general dynamical model for the energy exchange between the different scales in the system, 
as shown next. 

Recall that any large subsystem characterized by a temperature Tj is in contact with 
an even larger reservoir at temperature Since these temperatures differ, "heat" will 

flow between the two subsystems in accordance with Newton's law of cooling, so as to try 
to bring Tj close to Tj-\. In addition, there will be fluctuations in Tj of a random nature 
owing to the intermittency of the energy flow. Furthermore, the equations governing the 
temperature fluctuations must be invariant by a change, T — > AT, in temperature scale and 
ensure that the temperatures remain nonnegative. It then follows from these requirements 
that the temperature dynamics must be of the following form: 

dT 

= -H(T j -T j - 1 ) + g j (T J ,T j -.i)t i (t), (4) 

for j = 1, n, where the parameters /ij 1 correspond to the characteristic times of the prob- 
lem, the functions gj describe the noise amplitudes, and £j(t) denote mutually independent 
Gaussian white noises. Physically, the stochastic term in Eq. represents an effective 
coupling with the large-scale structures which accounts for intermittency The specific 
form of the function gj may depend on the system considered (see below), but it must pos- 
sess the following general properties: i) gj(XTj, A7}_i,) = \gj(Tj,Tj_i), on account of the 
invariance under change of temperature scale, and ii) gj(0, = 0, so as to ensure that 
the temperatures remain positive for all times (if they are initially positive). To see this, 
notice that if Tj = at some time, then Eq. implies that dTj/dt = fijTj_x > and so Tj 
never becomes negative. It is important to note that, irrespective of the form of gj, Eq. (T4J) 



implies that the internal reservoirs all have the same average temperature in the stationary 
regime, i.e., lim t ^ 00 (T J ) = T , for all j, as can be readily verified. 

We shall assume here that g 3 is a linear function of Tj, in which case Eq. (J3J) becomes 

dT 

~g = ~H(Tj - + VjT&Q), ( 5 ) 

where a 3 is a positive constant. This choice is not as restrictive as it seems, describing 
a rather general class of systems, as indicated below. Next we make use of the separa- 
tion of time scales, i.e., /x' 1 <C V>n-i "C • • • <S fJ-i 1 , in Eq. ([5]) to obtain the equilibrium 
distribution, f(T n ), of temperature in the innermost reservoir. In light of the time scale 
separation, it is safe to assume that over the characteristic time that it takes for the tem- 
perature T n _i of the surrounding reservoir to change appreciably, the temperature T n will 
relax to a quasi-stationary regime described by a conditional distribution, /(T„|T n „ 1 ), ob- 
tained from Eq. ([5]) for j = n with T n _i fixed. The marginal distribution for T n can 
then be written as a superposition of distributions /(T n |T n _i) with different values of T n _i: 
f(T n ) = J °° /(T n |T n _ 1 )/(T„_ 1 )<iT n „ 1 . Implementing this procedure recursively up to the 
outermost internal reservoir, one obtains 



poo poo n 

f^= / ••■ / Uf(T j \T j - 1 )dT 1 ---dT n -i. 
Jo Jo j=l 



(6) 



The distribution f(T 3 



Tj _i ) , for a given j, can easily be found by solving the stationary 



Fokker-Planck equation 16L Il7j associated with Eq. ([5]), holding Tj_i fixed. This yields an 



inverse gamma distribution 
where 

«i = (8) 
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If the system displays scale invariance one has aj = a, so that the distribution /(I}|7}_x) is 
identical across scales. That the parameter a 3 may become independent of scale is physically 
reasonable given that both \i 3 - and a 3 increase with j. (The latter follows from the fact 
that intermittency is stronger at smaller scales.) For scale invariant systems these two 
dependencies cancel out. 

With f(T n ) thus determined, let us now return to Eq. ([3]). To make further progress one 
needs to know the large-scale partition function Z (T) in order to carry out the integration 



over the variable T n . Let us consider the rather general case where the density of states, 
g{E), is a homogeneous function, that is, 

g{E)<xEP-\ 7 >0, (9) 

which implies that 

Z (T) oc {kT)\ (10) 

This relation describes several important classes of systems, such as: i) non-relativistic ideal 
gases, where 7 = //2, with / being the number of degrees of freedom of the system; and 
ii) systems where the energy is quadratic in the momenta and coordinates, in which case 

7 = /- 

Substituting Eq. (TTOj) into Eq. fl3J yields 

P r (ei) oc [°° (-^rYexp f(T n )dT n . (11) 



jo \ kT n J \ kT n j 
After inserting Eqs. (jSJ) and (J7J) into Eq. (jlll) . and performing a sequence of changes of 
variables of the type x n = a n T n _i/T n , one can show that the resulting multidimensional 
integral can be expressed in terms of known higher transcendental functions: 

PM) = ^-^tt n F (ai + 7 + 1, a n + 7 + 1; -fc), (12) 



where n F (ai, ...,a n ; —z) is the generalized hypergeometric function of order (n, 0) [18j . 
whose integral representation is given by 

n F (ai, ...,a n ) -z) = 

POO POO 

/ •■■/ e- x ^ z dX ai ( Xl )---dX an (x n ), (13) 
Jo Jo 

with d\ a (x) denoting the so-called Euler measure [19| : 

d\Jx) = 4re" I i tt " 1 ^. (14) 
r(a) 

1 1<=1 



In Eq. fll2]) the parameter j3 n is 
and the small-scale partition function Z n is given by 
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The generalized hypergeometric (GHG) distribution given in Eq. ( !T2|) has several interest- 
ing properties that are worth summarizing here. First, note that (E) r = / °° EP r (E)g(E)dE = 
7fcT . (This relation follows from the fact that {E) r = 7/c(T ra ) and (T n ) = T .) Thus, the 
equipartition theorem continues to hold on all scales in our multicanonical formalism. In this 
sense, the subsystems at different scales can be said to be in thermal equilibrium with one 
another at temperature To, even though the distribution law at small scales (i.e., for n > 1) 
differs from the usual BG distribution. Higher moments of the GHG distribution can also be 
readily computed: (E p ) r = (kT ) p r< ^^ YYi=i Yl^=i (^t]j- O ne t nen sees that, contrary to 
the canonical case, the multicanonical distribution P n (E) is not completely specified by its 
mean, which determines only the temperature To. Knowledge of the higher moments is nec- 
essary to determine the parameters cti. If the system displays scale invariance, i.e., «j = a, 
the value of a is determined by the second moment: (E 2 ) r = 7(7 + l)(/cT ) 2 [a/ (a — 1)]". 

Another important property of the GHG distribution is that it exhibits power-law tails 
of the form: P r {s) oc s^( a+ "/ +1 \ for e — > 00, as follows from the asymptotic expansion 



201 ] of the function n F (ai, ...,a n ; —x) for «j = a. It is also worth pointing out that the 
first two members of the family n F yield elementary functions, namely, o-^o (x) = exp(x) 
and iT (l/(g— l),x) = exp g (x/(q — 1)), where exp q (x) is the g-exponential: exp ? (x) = 
[1 + (1 - q)x] 1/(1 - q) . The GHG distribution with n = thus recovers the BG distribution, 



whereas for n = 1 it gives the g-exponential or Tsallis distribution {3]. For complex systems 
with more than two characteristic time scales GHG distributions of higher order are required. 

The multiscale formalism presented above can be readily extended to describe (statisti- 
cally stationary) fluctuations in highly-driven dissipative systems, such as fully-developed 
turbulence Although turbulent flows are out-of-equilibrium systems, the small-scale 
turbulence at high Reynolds numbers can be described in terms of an equilibrium theory, 



as first pointed out by Kolmogorov 



211 ] . This means that the small eddies in the range 



re local 
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23]. 



r<L, where L is the integral scale at which energy is injected, quickly adjust to t 
conditions of the mean flow and are therefore in approximate statistical equilibrium 
Furthermore, in the inertial subrange (i.e., for rj <^ r <^ L, where n is the Kolmogorov scale 
at which viscous effects become relevant), energy is transferred from large eddies to smaller 
ones with essentially no dissipation. In our formalism, the local-equilibrium condition is 
contained in Eqs. ([3]) and (|6]), whereas energy conservation corresponds to the fact that 
(Tj) = T in the stationary regime. 



Now we wish to apply our multiscale formalism to Lagrangian turbulence, where one 
is concerned with the dynamical properties of individual fluid particles. In Lagrangian 
turbulence, intermittency manifests itself as a change in shape of the PDF of velocity time 
increments with the time lag. Let us then consider time increments, 5 T v — v(t + r) — v (i), 
of one component v of the Lagrangian velocity. Here we take r = Tl/2 u , where Ti is the 
integral time scale which is related to the large-eddy turnover time. The fluctuations in the 



velocity increments can be modelled 

d(5 T v 



10| by a Langevin equation of the type 

-rti T v + ££(*), (17) 



dt 

where the "friction coefficient" T is assumed to be constant, but the noise amplitude S is 
allowed to fluctuate in a slow time scale as compared to the relaxation time T -1 , which 
is of the same order of magnitude as the Kolmogorov time r rj . It then follows that over 
short time scales the velocity fluctuations reach a quasi-equilibrium described by a Gaussian 
distribution, 

»^^f exp (-W)' (18) 

with variance cr 2 . = S 2 /T, which is assumed to be proportional to the fluctuating energy 



dissipation rate, e, times the time leg r 10|, |24j. In the context of our multicanonical 
formalism, <j 2 T plays the role of the fluctuating temperature T n , compare Eqs. (ED) and (TTBI) . 
so that the marginal distribution P(5 T v) of velocity time increments can be written in a 
form equivalent to Eq. (fT2l . only replacing e with (8 T v ) 2 and setting 7 = 1/2 (corresponding 
to one degree of freedom, since only one velocity component of the Lagrangian particle is 
considered). One then obtains that P(S T v), normalized to unit variance, is given by the 
following GHG distribution: 



P(S T v) 



r(q + 3/2) 
_a 1 / 2 T(a + 1) 



t F (a + 3/2, a + 3/2; -^p), (19) 



where we have set oti = a. An earlier derivation of the distribution (|T9|) was given in Ref. 25 1 



in the context of Eulerian turbulence, where starting from the scale-by-scale energy budget 



equation \l\ we proposed a set of SDE's similar to Eq. (J5]). As noted above, this distribution 
is but a particular case of the more general GHG distribution given in Eq. ( |T2l) . 

An application of the GHG distribution ( TT9|) to Lagrangian turbulence is shown in Fig. [TJ 
In this figure the circles represent the PDF of acceleration measurements on a turbulent 
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FIG. 1. Distribution of accelerations on turbulent water flow at R\ = 690 (circles) and theoretical 
curve (solid line) for a GHG distribution with N = 7 and a = 2.67. 



water flow (R\ = 690) performed in Prof. Bodenschatz's group |4|; for details about the 



experiments see Ref. |26|. Because the acceleration, a, was computed from the position 
measurements by a filtering procedure of width ~ t v {4], it is safe to assume that a is 
proportional to 5 T v for r m t v . The number n of scales can then be estimated as: n = 
log 2 (T L /r v ) = log 2 (R\/ \/lE)- Thus, for R x = 690 one finds n = 7. The parameter a, on the 
other hand, can be estimated by matching the fourth moment of the empirical distribution, 
which yields a = 2.67. Superimposed with the experimental data in Fig. [I] is the plot of the 
GHG distribution for n = 7 and a = 2.67. The agreement between the theoretical curve 
and the data is remarkable. The dependence of a on R\ can be obtained by noting that 
the acceleration flatness, F = (a 4 ) /(a 2 ) 2 = 3[a/(a — l)] n , initially increases with R\ but 
then seems to level off for R\ > 500 [26]. If this tendency holds for R\ — > 00, one must 
have a = C^ 1 log 2 R\, where C is a constant, so that a oc n, which then yields F = 3e c as 
n — t- 00. In this limit, the GHG distribution recovers the log-normal model widely used in 



turbulence 



251 ] . where it has been conjectured [27j that C is a universal constant (C ~ 3) for 



R\ — > 00. Our results show however that the log-normal model holds only asymptotically 
as R\ — > 00, whereas for finite R\ the GHG distribution should apply. 

As a concluding remark, we note that the distribution ([7]) can be derived from a maximum 
entropy principle by extending the arguments used in Refs. 28|, |29[ for the case n — 1. More 
details will be published elsewhere; here it suffices to say that in this approach the parameter 
a appears as a Lagrange multiplier and the connection between Eqs. ([7]) and (Q becomes 



more apparent. It is perhaps also worth noting that an alternative derivation of Eq. (J7J) 
can be given on the basis of Bayesian inference [30] . Other multiscale systems are currently 
under investigation. 
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